Method for analyzing nonlinear restoring force characteristic with hysteresis of machine structure system

ABSTRACT

The invention discloses a method for analyzing nonlinear restoring force characteristic with hysteresis of a machine structure system, capable of easily and high-accurately performing analysis of nonlinear vibration characteristic of the machine structure system having a displacement strain characteristic dependent on a past history. When analysis is performed for the nonlinear restoring force characteristic of the machine structure system, having force-displacement hysteresis, approximate mode analysis is carried out by representing the nonlinear restoring force characteristic having the history by an equivalent linear characteristic.

BACKGROUND OF THE INVENTION

1. Field of the Invention

This invention relates to an approximate mode analysis method (method ofapproximate modal analysis) to analyze nonlinear restoring forcecharacteristic with force-displacement hysteresis of a machine structuresystem. More particularly, the invention relates to the method ofanalysis for improving the vibration characteristic of the machinestructure system by reconstructing a hysteretic characteristic using thedata of a discrete load strain characteristic obtained experimentally bya vibration test or the like, alternatively, simulatively by a finiteelement method (FEM) analysis or the like, and by supplying an outputaccording to the obtained data with respect to an arbitrary input.

2. Description of the Related Art

The mode analysis method has been established for its technology as amethod for analyzing the vibration of the machine structure system, andhas also been made commercially available as software, and frequentlyused. However, while the mode analysis method itself is based on atheory assuming a linear system, an actual machine structure systemnormally includes a plurality of nonlinear characteristics, in a boltjoined part, a sliding surface part of a baring or the like, arotational part, a rubber pad of an engine mount as shown in FIG. 19, ajoint of a robot, and so on.

With regard to a case including a nonlinear characteristic having nohysteretic characteristic, presented is a method for approximatelyachieving mode analysis by representing the nonlinear characteristic byan equivalent linearization method. This method has contributed to therigidity designing and the working accuracy improvement of a machinetool, the support characteristic evaluation of a piping system andreflection in an aseismic designing method, optimal designing in supportat the engine mount, and so on.

However, a restoring force characteristic typically represented by therubber pad exhibits hysteretic nonlinearity. Conventionally, it wasdifficult to reduce an error when the estimated value obtained byestimating the frequency response characteristic of the machinestructure system including such a nonlinear characteristic in adesigning stage was compared with the frequency response characteristicobtained by an experiment.

When a vibration characteristic is requested by considering a structureof element, a characteristic for connecting elements together, and soon, in the designing stage, generally, analysis is carried out by usingthe finite element method. However, for a vibration characteristic,analysis is carried out mainly based on time history, and to obtain afrequency response characteristic, time history response must berepeated for each frequency step. Thus, a great deal of computing timeis required even for obtaining the response of the linear system Inaddition, it is also possible to obtain a frequency responsecharacteristic from the equation of motion. Regarding the case includingthe nonlinear characteristic, however, no general and simple methodshave been presented to obtain the vibration characteristic of the entiremachine structure system by presenting the nonlinear characteristic.

By carrying out the vibration experiment of the machine structuresystem, a frequency response characteristic including the abovecharacteristic is obtained. In the mode analysis method, this is calledan experimental mode analysis. According to this method, a frequencyresponse characteristic is obtained by an experiment, and the systemcharacteristic is improved by estimating an intrinsic frequency, adamping constant, a mass, spring, damping coefficients, and so on.

However, if the nonlinear characteristic is included, a frequencyresponse characteristic is obtained in a biased manner from a linearcharacteristic. Consequently, because of a biased frequency, distortionof a characteristic in the vicinity of the intrinsic frequency, and soon, it is difficult to estimate an intrinsic frequency, and generally, adamping constant is excessively evaluated. As a result, wrong treatmenthas been taken for an improvement of the characteristic, or no propertreatment has been taken therefor.

With regard to nonlinear restoring force characteristic withouthysteresis, a plural nonlinear simultaneous equation can be constructedby applying the equivalent linearization method to amulti-degree-of-freedom system equation of motion. It has beenconsidered possible to find a solution of the equation, i.e., thefrequency response characteristic of the multi-degree-of-freedom systemincluding a plurality of nonlinear characteristics. In addition, byusing a building block method (BB method) connected with the finiteelement method to connect the vibration characteristic of the devicestructure by the nonlinear restoring force characteristic, it has beenconsidered possible to approximately carry out the mode analysis of themulti-degree-of-freedom nonlinear vibration system on the extension ofthe conventional mode analysis method.

However, the foregoing method was unable to deal with the systemincluding the nonlinear restoring force characteristic with hysteresis,as typically seen in the rubber pad used when an automotive engine wasmounted on a car body. In addition, the nonlinear characteristic of abearing and the like was often represented by hysteretic nonlinearcharacteristic, and it was impossible to deal with this characteristic.

In the vibration characteristic analysis of the mechanical system, forthe method of representing a hysteretic characteristic in a rigid part,a model for representing the hysteretic characteristic by a combinationof divisional straight lines, such as bilinear model, a trilinear modelor the like, a hysteretic characteristic model for representing it by anumerical expression, such as Ramberg-Osgood type, and so on, have beenused. However, considering the simulation of the nonlinearcharacteristic of an actual machine based on the characteristics of suchmodels, problems described below inevitably occur.

Each of the bilinear and trilinear models is a method for representing ahysteretic characteristic by the combination of straight lines. In thesemethods, the hysteretic characteristics to be represented are consideredto be two and three inclined straight lines in sections respectively inthe bilinear and trilinear models, and the shapes of the hystereticcharacteristics are simulated by joining the straight lines oversections. A problem inherent in each of these methods is that since thehysteretic characteristic to be presented takes a complex shape, itcannot be sufficiently represented by the sectional combination of twoto three straight lines.

In addition, since the hysteretic characteristic to be represented isconsidered to be two to three inclined straight lines, the inclinationis changed piecewise. Consequently, even when the hystereticcharacteristic to be represented takes a smooth shape, the sectionalchange of the inclination causes an output from the model to becomeunsmooth with respect to an arbitrary input.

There is also a problem inherent in the case of representing thehysteretic characteristic by the bilinear model. That is, when sectionalstraight line regression is made with respect to the hystereticcharacteristic to be represented, depending on the position ofinclination changing to be set, an obtained result varies even if asimilar hysteretic characteristic is represented.

On the other hand, the hysteretic characteristic model represented bythe Ramberg-Osgood type is a method for representing a hystereticcharacteristic by a polynomial. This method is designed to simulate theshape of the hysteretic characteristic by deciding maximum and minimumcoordinates based on a skeleton curve, and by connecting upward anddownward curves from the coordinates. To represent the hystereticcharacteristic by this method, it is necessary to match the hystereticcharacteristic to be represented with the hysteretic characteristicmodel by using a least square method or-the like. However, when thehysteretic characteristic obtained from the vibration test, FEManalysis, or the like has a complex shape, it is difficult to set aparameter to match the characteristic in detail with the curve-matchedhysteretic characteristic model.

Moreover, in the system represented by the hysteretic characteristicmodel of the Ramberg-Osgood type, for performing approximate vibrationcharacteristic analysis by the equivalent linearization method, it isnecessary to calculate Fourier series of an output with respect to anarbitrary input, as an approximation of a nonlinear factor. However,depending on a parameter for representing this model, it may bedifficult to analytically calculate Fourier series.

A difference in equivalent rigidity and equivalent damping depending onthe presence of a hysteretic characteristic is that while there isneither phase delay nor advancement in a relation between an input andan output in the case of absence of hysteretic characteristics, in thecase of the presence of a hysteretic characteristic, a phase relation isincluded in a relation between an input and an output, and this need beconsidered when united with the equation of motion.

With regard to the representation of nonlinear hysteretic characteristicby a describing function (equivalent linearization method in the fieldof automatic control), in the field of an electric system, when thecharacteristics of a current and a magnetic field were taken intoconsideration, it was considered possible to numerically represent thecharacteristics by the describing function based on actualcharacteristics. However, for the application to themulti-degree-of-freedom system and the representation of the frequencyresponse characteristic as mode analysis as seen in the machinestructure system, no particular solution of a plural nonlinearsimultaneous equation having a phase characteristic from the beginninghas been presented, except for an area targeted from a mathematicalinterest.

For performing vibration characteristic analysis of the nonlinear systemhaving a displacement strain characteristic dependent on a past history,it was necessary to apply the straight line regression, the least squaremethod, or the like to the foregoing hysteretic characteristic model.Thus, it was difficult to match, in detail, the shape of the hystereticcharacteristic model used for analysis with the data of the hystereticcharacteristic obtained from the vibration test, the FEM analysis, orthe like to be applied, and this point was a problem when vibrationcharacteristic analysis was carried out.

SUMMARY OF THE INVENTION

The present invention was made to advantageously solve the foregoingproblems. In accordance with the invention, a method of analysis isprovided for the representation of the vibration characteristic of amachine structure system having nonlinear restoring force characteristicwith hysteresis by a frequency response characteristic, the method ofanalysis making approximate mode analysis enable to be performed byrepresenting the vibration characteristic numerically by equivalentrigidity and equivalent damping based on the nonlinear characteristic ofan actual system, and by uniting such with a multi-degree-of-freedomequation of motion and obtaining a solution of the equation. This methodcontributes to the evaluation and improvement of the vibrationcharacteristic of the machine structure system. As a result, it ispossible not only to predict the characteristic of the system in adesigning stage, but also to estimate a vibration system parameter moreaccurately and rationally from the vibration characteristic obtained bythe vibration experiment of the actual system.

According to the invention, without using the method of representing ahysteretic characteristic based on the application of hystereticcharacteristic data obtained from a vibration test or the like to acertain model, an output is made according to the obtained hystereticcharacteristic with respect to an arbitrary input, based on the obtainedhysteretic characteristic data. For such a purpose, Preisach modeldesigned for recording a magnetic characteristic in an electric andelectronic field is used to make a hysteretic characteristic modelregarding the load strain of the mechanical system. Moreover, amechanical vibration characteristic analysis including a testcharacteristic is performed by using this hysteretic characteristicrepresenting method.

When a hysteretic characteristic model obtained from the vibration testor the like is made by the Preisach model, since this model includesmany relay elements, it is necessary to measure a minor loopconstituting this hysteretic characteristic model. However, themeasurement of the minor loop of the hysteretic characteristic by thevibration test may be difficult. In addition, there is a problem whenthe hysteretic characteristic model obtained by the invention is usedfor vibration characteristic analysis. That is, when an equivalentlinearization method is used because the hysteretic characteristic isnumerically represented, it is impossible to analytically obtainequivalent rigidity obtained by Fourier transformation. Thus, accordingto the invention, equivalent rigidity is numerically obtained, and usedfor analysis.

According to the invention, when a solution is found to the nonlinearequation of motion of the system by Newton-Raphson method, data isobtained on the equivalent rigidity of a nonlinear factor, and this datais interpolated by a piecewise polynomial such as a spline function orthe like. The spline function represented by the piecewise polynomialcan be differentiated. Thus, Jacobian matrix is formed whenNewton-Raphson method is applied. This process develops an output to thehysteretic characteristic in Fourier series by using a solution obtainedfor each repetition, making it possible to achieve high efficiency forthe method of obtaining equivalent rigidity and equivalent damping.

According to the invention, the vibration characteristic analysisincluding the hysteretic characteristic data obtained from the vibrationtest, FEM analysis, or the like, can be performed. Thus, it is possibleto perform analysis more closely related to an experiment. Moreover,according to the hysteretic characteristic analysis method of theinvention, even a nonlinear characteristic having no history can beanalyzed from data obtained from the vibration test, the FEM analysis,or the like by a similar method.

The foregoing arrangements enable a frequency response characteristic tobe clarified by taking, into the equation of motion, a nonlinearcharacteristic existing in a connected part between element structuresin a machine tool or the like, a nonlinear characteristic provided in abearing, a sliding surface part or the like, a nonlinear characteristicof a support component used for engine suspending, a nonlinearcharacteristic generated in a robot joint, and so on. Accordingly, theimprovement of the characteristic and dealing with the characteristicbased on the above considerations can be rationally pursued. Comparedwith the approximate mode analysis based on a consideration given to thenonlinear characteristic having no history, more accurate estimation canbe made quickly by numerically representing the hystereticcharacteristic by equivalent rigidity and equivalent damping andconsidering the same based on an actual characteristic.

Therefore, in an age where competition is severe in an automobileindustry or the like, among nations, and between corporate groups, bysolving the existing problems, it is possible to realize the promptstarting of production, and to quickly deal with problems which arise.

BRIEF DESCRIPTION OF THE DRAWINGS

FIG. 1 is a flowchart showing a process of a method of analysisaccording to an embodiment of the invention;

FIGS. 2(a) to 2(c) are views, each thereof illustrating a configurationof a system of analysis used for the method of analysis of theembodiment;

FIG. 3 is a characteristic view showing an example of nonlinearrestoring force characteristic with hysteresis;

FIGS. 4(a) to 4(c) are illustrations: FIG. 4(a) showing a minimum factorof a hysteretic characteristic represented by a discrete value; FIG.4(b) showing a hysteretic characteristic represented by a collection ofminimum factors; and FIG. 4(c) showing an area, in which a distributionfunction is present on an x_(d)-x_(i) plane;

FIG. 5 is a view illustrating Preisach distribution function representedby a discrete value seen in a three-dimensional space;

FIGS. 6(a) to 6(c) are views showing a change in a load when adisplacement increases;

FIGS. 7(a) to 7(c) are views showing a change in a load when adisplacement decreases;

FIGS. 8(a) and 8(b) are views showing an area of matched addition of adistribution function table by discrete values;

FIG. 9 is a flowchart showing a process for representing a hystereticcharacteristic by a discrete value;

FIG. 10 is a flowchart showing an algorithm for making a distributionfunction table from vibration test data;

FIG. 11 is a view illustrating a data acquisition point when making thedistribution function;

FIG. 12 is a view illustrating a method for making a hystereticcharacteristic model (amplitude X₁) by processing the vibration testdata;

FIGS. 13(a) and 13(b) are views illustrating a method for making ahysteretic characteristic model (amplitude X₂) by processing thevibration test data;

FIG. 14 is a view illustrating Preisach distribution function formedfrom numerical model data, seen in a three-dimensional space;

FIG. 15 is a view illustrating hysteretic characteristic data formed bya numerical model, and a hysteretic characteristic by Preisach model;

FIG. 16 is a flowchart showing a process for obtaining a solution of amulti-degree-of-freedom nonlinear simultaneous equation;

FIG. 17 is a characteristic line view showing an example of anequivalent rigidity coefficient k_(e)(X) and a differential valuek_(e)′(X) thereof;

FIG. 18 is a characteristic line view showing an example of a frequencyresponse characteristic of single-degree-of-freedom system withhysteresis; and

FIG. 19 is a view illustrating an engine shake model having an enginemount as an example of a machine structure system with nonlinearrestoring force characteristic (K_(NL)).

DESCRIPTION OF THE PREFERRED EMBODIMENT

Next, the preferred embodiment of the present invention will bedescribed with reference to the accompanying drawings.

FIG. 1 is a flowchart showing a procedure of an analysis methodaccording to an embodiment of the present invention; and each of FIGS.2(a) to 2(c) is a view exemplifying a configuration of a system ofanalysis used for the analysis method of the embodiment. To give anoutline of the analysis method of the embodiment, first, as shown inFIG. 2(a), data is obtained on a restoring force characteristic withhysteresis of a load output F with respect to an input displacement X bymaking an excitation experiment for a test piece 2 such as a rubber padto be analyzed at the exciter 1 of an excitation testing device, and/orsubjecting the test piece to simulation analysis (simulation by acomputer) such as FEM analysis, and so on. The obtained data is thenentered through an oscilloscope 3 to a personal computer (PC) 4. Then,as shown in FIG. 1, Preisach model is made by data processing of the PC4, equivalent rigidity and an equivalent damping coefficient areobtained by an equivalent linearization method, and then the equivalentrigidity and the equivalent damping coefficient are substituted for anequation of motion of multi-degree-of-freedom system. By solvingnonlinear simultaneous equations, for example, a frequency responsecharacteristic like that shown in FIGS. 2(b) and 2(c) is obtained.

Specifically, according to the analysis method of the embodiment, when ahysteretic characteristic model is constructed, the data of thehysteretic characteristic regarding a load and strain obtained from avibration test, FEM analysis or the like is processed by being dividedinto a case where the input displacement is increased and a case wherethe input displacement is decreased. In other words, the displacement Xis divided into a displacement x_(i) at the time of increase and x_(d)at the time of decrease. A relation between the displacement X and aload F is represented based on these two variables of displacementsx_(i) and x_(d), and a hysteretic characteristic model is made asdiscrete values by obtaining a change ΔF in the load with respect tochanges in the increase and the decrease of the displacement.

As shown in FIG. 3, regarding to the relation between the displacementand the load, if the load is represented by the following continuousfunction when the displacement is increased,

F=j(x)  (1)

then the increase of the load F with respect to the displacement x isrepresented as follows: $\begin{matrix}{\frac{F}{x} = {j^{\prime}(x)}} & (2)\end{matrix}$

The equation (2) represents a slope of the relation between thedisplacement and the load. Similarly, when the displacement isdecreased, if the relation between the load and the displacement isrepresented by the following continuous function,

 F=k(x)  (3)

then the reduction of the load F with respect to the displacement x isrepresented as follows: $\begin{matrix}{\frac{F}{x} = {k^{\prime}(x)}} & (4)\end{matrix}$

Here, the displacement x is divided into displacement at the time ofincrease and displacement at the time of decrease respectively set asindependent variables, and the displacements at the time of increase andat the time of decrease are respectively denoted by x_(i) and x_(d). Achange in the load F with respect to the changes of the displacementsx_(i) and x_(d) can be represented as follows: $\begin{matrix}{\frac{\partial^{2}F}{{\partial x_{i}}{\partial x_{d}}} = {\eta \left( {x_{i},x_{d}} \right)}} & (5)\end{matrix}$

By integrating the equation (5) with respect to x_(d), a result is theslope of the load with respect to the displacement at the time ofdisplacement increase of the equation (2). Similarly, by integrating theequation (5) with respect to x_(i), a result is the slope of the loadwith respect to the displacement at the time of displacement decrease ofthe equation (4). By integrating these results in the range of theincrease and decrease changes of the displacement, then a load withrespect to a given displacement is outputted. This η(x_(i), x_(d)) iscalled as Preisach distribution function.

When considered by discrete value, the Preisach distribution functionbecomes as follows. Consideration is given to a characteristic like thatshown in FIG. 4(a), represented by providing three values, i.e.,displacement values x_(i), x_(d) and a load value f. This is representedby F=f(x_(i), x_(d)).

(1) When a displacement satisfies X<x_(i), F=−f/2 is established.

(2) When the displacement is further increased, F=f/2 is established atthe point of X=x_(i), and this state is maintained.

(3) When the displacement is decreased, F=f/2 is established withX<x_(d).

(4) When the displacement is further decreased, F=−f/2 is established atthe point of X=x_(d), and this state is maintained.

Now, assuming that one having the foregoing characteristic is a minimumfactor of a hysteretic characteristic, and as shown in FIG. 4(b), thehysteretic characteristic is considered to be a collection of theminimum factors, represented as a square type.

Regarding the above hysteretic characteristic, considering that noreverse movements occur in routes during the increase and the decrease,if the maximum and minimum values of the displacement x of thecharacteristic are X_(s) and −X_(s), then the defined domain of theabove η(x_(i), x_(d)) becomes as follows:

x _(i) ≦x _(d) , |x _(i) |≦X _(s) , |x _(d) |≦X _(s)  (6)

When the Preisach distribution function is seen on a plane x_(i)-x_(d),the defined domain of this distribution function becomes a triangulararea like that shown in FIG. 4(c). FIG. 5 shows the Preisachdistribution function seen in a three-dimensional space. (Change in loadat time of displacement increase)

Now, the distribution function will be described by using the planex_(i)-x_(d). When the displacement is increased from X to X+ΔX(ΔX>0),the sum of all the distribution functions η(x_(i), x_(d)) when thedisplacement x_(i) at the time of increase is between X and X+ΔX becomesan increase ΔF of the load, which is represented as follows.$\begin{matrix}{{\Delta \quad F} = {\int_{- X_{s}}^{X_{s}}\quad {{x_{d}}{\int_{X}^{X + {\Delta \quad X}}{{\eta \left( {x_{i},x_{d}} \right)}\quad {x_{i}}}}}}} & (7)\end{matrix}$

FIGS. 6(a) to 6(c) show a change in the integration domain of thedistribution function at the time of displacement increase regarding thehysteretic characteristic of amplitude X_(s). At the position of (1) inFIG. 6(a), the load F takes a minimum value F_(min), which isrepresented as follows. $\begin{matrix}{F_{\min} = {{- \frac{1}{2}}{\int_{- X_{s}}^{X_{s}}\quad {{x_{d}}{\int_{- X_{s}}^{X_{s}}{{\eta \left( {x_{i},x_{d}} \right)}\quad {x_{i}}}}}}}} & (8)\end{matrix}$

The sum of the F_(min) and ΔF obtained by integrating η(x_(i), x_(d)) inthe displacement increase region becomes a load F. $\begin{matrix}{F = {F_{\min} + {\int_{- X_{s}}^{X_{s}}\quad {{x_{d}}{\int_{X}^{X + {\Delta \quad X}}{{\eta \left( {x_{i},x_{d}} \right)}\quad {x_{i}}}}}}}} & (9)\end{matrix}$

A hatched part for integration is increased from the position of −X_(s)of (1) to 0 of (2) in FIG. 6(b), and increased to X_(s) of (3) in FIG.6(c), becoming a maximum value F_(max) of F. (Change in load at time ofdisplacement decrease)

When the displacement is decreased from X to X−ΔX_(s) the sum of all thedistribution functions η(x_(i), x_(d)) when the displacement x_(d) atthe time of decrease is between X and X−ΔX becomes a load reduction ΔF,which is represented as follows. $\begin{matrix}{{\Delta \quad F} = {{\int_{X}^{X - {\Delta \quad X}}\quad {{x_{d}}{\int_{- X_{s}}^{X_{s}}{{\eta \left( {x_{i},x_{d}} \right)}\quad {x_{i}}}}}} = {- {\int_{X - {\Delta \quad X}}^{X}\quad {{x_{d}}{\int_{- X_{s}}^{X_{s}}{{\eta \left( {x_{i},x_{d}} \right)}\quad {x_{i}}}}}}}}} & (10)\end{matrix}$

FIGS. 7(a) to 7(c) show a change in the integration domain of thedistribution function at the time of displacement decrease. The load Ftakes a maximum value F_(max) at the position of (4) in FIG. 7(a), whichis represented as follows. $\begin{matrix}{F_{\max} = {\frac{1}{2}{\int_{- X_{s}}^{X_{s}}\quad {{x_{d}}{\int_{- X_{s}}^{X_{s}}{{\eta \left( {x_{i},x_{d}} \right)}\quad {x_{i}}}}}}}} & (11)\end{matrix}$

The sum of the F_(max) and ΔF obtained by integrating η(x_(i), x_(d)) ofthe displacement decrease region becomes a load F. $\begin{matrix}{\quad {F = {{F_{\max} + {\int_{X}^{X - {\Delta \quad X}}\quad {{x_{d}}{\int_{- X_{s}}^{X_{s}}{{\eta \left( {x_{i},x_{d}} \right)}\quad {x_{i}}}}}}} = {F_{\max} - {\int_{X - {\Delta \quad X}}^{X}\quad {{x_{d}}{\int_{- X_{s}}^{X_{s}}{{\eta \left( {x_{i},x_{d}} \right)}\quad {x_{i}}}}}}}}}} & (12)\end{matrix}$

A hatched part is reduced from the position of X_(s) of (4) in FIG. 7(b)to 0 of (5), and then reduced to −X_(s) of (6) in FIG. 7(c), becoming aminimum value F_(min) of F.

The foregoing is a concept of the hysteretic characteristicrepresentation by the Preisach model.

The representation of the foregoing by discrete values is as follows.That is, it is assumed that the number of divisions in the distributionfunction represented by discrete values is M, and amplitude X as aninput maximum value is present in the Nth division. When thedisplacement is increased by ΔX from X at the Nth division, if theposition of the division having ΔX is N_(x), a change in the load withrespect to the increase of the displacement may be like that shown inFIG. 8(a), and represented as follows. $\begin{matrix}{F_{\min} = {{- \frac{1}{2}}{\sum\limits_{x_{i} = 0}^{M - 2}\quad {\sum\limits_{x_{d} = 0}^{M - {2N}}\quad {\eta \left( {x_{i},x_{d}} \right)}}}}} & (13) \\{F = {F_{\min} + {\sum\limits_{x_{i} = N}^{N_{x} - N}\quad {\sum\limits_{x_{d} = N}^{M - {2N}}\quad {\eta \left( {x_{i},x_{d}} \right)}}}}} & (14)\end{matrix}$

Similarly, when the displacement is decreased by ΔX from X at the Nthdivision, if the position of the division having ΔX is Nx, a change inthe load with respect to the decrease of the displacement may be likethat shown in FIG. 8(b), and represented as follows. $\begin{matrix}{F_{\max} = {\frac{1}{2}{\sum\limits_{x_{i} = 0}^{M - {2N}}\quad {\sum\limits_{x_{d} = 0}^{M - {2N}}\quad {\eta \left( {x_{i},x_{d}} \right)}}}}} & (15) \\{F = {F_{\max} + {\sum\limits_{x_{i} = N}^{M - {2N}}\quad {\sum\limits_{x_{d} = N}^{N_{x} - N}\quad {\eta \left( {x_{i},x_{d}} \right)}}}}} & (16)\end{matrix}$

The flowchart of the entire process of the hysteretic characteristicrepresentation by such discrete values is shown in FIG. 9.

Now, the process of representing the hysteretic characteristic will bedescribed in detail with reference to a flowchart shown in FIG. 10.

(Step 1)

Consideration is given to a plane x_(i)-x_(d) regarding the displacementx_(i) at the time of increase and the displacement x_(d) at the time ofdecrease. By a vibration test, FEM analysis or the like, regarding aloop of the hysteretic characteristic of a target test piece such as arubber pad, as shown in FIG. 11, data is obtained, in which theamplitude of the loop is increased to X₁, X₂, . . . , X_(n) at equalintervals ΔX. Then, the values of loads are obtained with respect to thedisplacement values −X_(n), −X_(n−1), . . . , −X₁, 0, X₁, X₂, . . . ,X_(n) of the hysteretic characteristic data. If there are no load dataat the displacement values −X_(n), −X_(n−1), . . . , −X₁, 0, X₁, X₂, . .. , X_(n), the load data is obtained by interpolating data.

(Step 2)

The above hysteretic characteristic data is divided into load data atthe time of displacement increase and load data at the time ofdisplacement decrease. Regarding the displacement increase at amplitudeX_(k) (k=1, 2, . . . , n) and loads at the time of displacementincrease, if a load at the minimum value −X_(k) of the amplitude isF_(−k, 2k), and a load at a maximum value X_(k) is F_(k, 2k), then 2k+1pieces of data are represented as follows:

x _(i) ^((k)) =[−X _(k) , −X _(k−1) , . . . , −X ₁ , X ₀ , X ₁ , X ₂ , .. . , X _(k)](X ₀=0)  (17)

F _(i) ^((k)) =[F _(−k, 2k) , F _(k, 1) , F _(k, 2) , . . . , F_(k, 2k)]  (18)

Loads at the time of displacement decrease with respect to displacementdecrease are represented as follows:

X _(d) ^((k)) =[X _(k) , X _(k−1) , . . . , −X ₁ , X ₀ , −X ₁ , −X ₂ , .. . , −X _(k)](X ₀=0)  (19)

F _(d) ^((k)) =[F _(k, 2k) , F _(−k, 1) , F _(−k, 2) , . . . , F_(−k, 2k)]  (20)

(Step 3)

A change ΔF_(i) in load is obtained with respect to a displacementincrease. The change ΔFi is represented as follows: $\begin{matrix}\begin{matrix}{{\Delta \quad F_{i}^{(k)}} = \quad \left\lbrack {F_{{ik},1},F_{{ik},2},\ldots \quad,F_{{ik},{2k}}} \right\rbrack} \\{= \quad \left\lbrack {{F_{k,1} - F_{{- k},{2k}}},{F_{k,2} - F_{k,1}},{F_{k,3} -}} \right.} \\{\quad \left. {F_{k,2},\ldots \quad,{F_{k,{2k}} - F_{k,{{2k} - 1}}}} \right\rbrack}\end{matrix} & (21)\end{matrix}$

A change ΔF_(d) in load is obtained with respect to a displacementdecrease. The change ΔF_(d) is represented as follows: $\begin{matrix}\begin{matrix}{{\Delta \quad F_{d}^{(k)}} = \quad \left\lbrack {F_{{dk},1},F_{{dk},2},\ldots \quad,F_{{dk},{2k}}} \right\rbrack} \\{= \quad \left\lbrack {{F_{{- k},1} - F_{k,{2k}}},{F_{{- k},2} - F_{{- k},1}},{F_{{- k},3} -}} \right.} \\{\quad \left. {F_{{- k},2},\ldots \quad,{F_{{- k},{2k}} - F_{{- k},{{2k} - 1}}}} \right\rbrack}\end{matrix} & (22)\end{matrix}$

(Step 4)

By use of the above load changing quantities with respect to thedisplacement changes, η(x_(i), x_(d)) is obtained by processing ΔF_(i)and ΔF_(d) sequentially from k=1 to k=n as follows. There are threeη(x_(i), x_(d)) at k=1, and represented as follows: $\begin{matrix}\begin{matrix}{{\eta^{(1)}\left( {x_{i},x_{d}} \right)} = \quad \left\lbrack {{\eta \left( {X_{0},X_{1}} \right)},{\eta \left( {{- X_{1}},X_{1}} \right)},{\eta \left( {{- X_{1}},X_{0}} \right)}} \right\rbrack} \\{= \quad \left\lbrack {F_{i12},{F_{d11} - F_{i12}},F_{d12}} \right\rbrack} \\{= \quad \left\lbrack {F_{i12},{F_{d11} - {\eta \left( {X_{0},X_{1}} \right)}},F_{d12}} \right\rbrack}\end{matrix} & (23)\end{matrix}$

There are seven η(x_(i), x_(d)) at k=2, which are represented as followsby using F_(i12), F_(i11), F_(d12), and F_(d11) obtained at k=1:

η⁽²⁾(x _(i) , x _(d))=[η(X ₁ , X ₂), η(X ₀ , X ₂), η(−X₁ , X ₂), η(−X₂ ,X ₂), η(−X₂ , −X ₁), η(−X₂ , X ₀), η(−X₂ , X ₁)]

=[F _(i24) , F _(i23) −F _(i12) , F _(i22) −F _(i11) , F _(i21)−[η(X ₁ ,X ₂)+η(X ₀ , X ₂)+η(−X ₁ , X ₂), F _(d22) −F _(d11) , F _(d23) −F _(d12), F _(d24)]  (24)

Similarly, there are 4n−1 pieces of η(x_(i), x_(d)) at k=n, which arerepresented as follows by using F_(in−1, 2n−2), F_(in−1, 2n−3),F_(in−1, 1), F_(dn−1, 1), F_(dn−1, 2n−3), F_(dn−1, 2n−2) obtained atk=n−1. $\begin{matrix}\left. \begin{matrix}{{{\eta \left( {X_{n - 1},X_{n}} \right)} = F_{{i\quad n},{2n},}}\quad} \\{{{\eta \left( {{- X_{n}},{- X_{n - 1}}} \right)} = F_{{dn},{2n},}}\quad} \\{{\eta \left( {{- X_{n}},X_{n - 1}} \right)} = {F_{i\quad {n1}} - \left\{ {{\eta \left( {X_{n - 1},X_{n}} \right)} + {\eta \left( {X_{n - 2},X_{n}} \right)} + \ldots + {\eta \left( {{- X_{n}},X_{n}} \right)}} \right\}}}\end{matrix} \right\} & (25)\end{matrix}$

 η^((n))(x _(i) , x _(d))=[η(X _(n−1) , X _(n)), η(X _(n−2) , X _(n)), .. . , η(−X _(n) , X _(n)), η(−X _(n) , X _(n−1)), . . . , η(−X _(n) , X_(n−1))]

=[F _(in, 2n) ,

F_(in, 2n−1) −F _(in−1, 2n−2) ,

F_(in, 2n−2) −F _(in−1, 2n−3) , . . . ,

F_(in2) −F _(in−1, 1) ,

F_(in1)−{(η(X _(n−1) , X _(n))+η(

X_(n−2) , X _(n))+. . . +η(−

X_(n) , X _(n))}, F _(dn, 2) −F _(dn−1, 1) , . . . ,

F_(dn, 2n−2) −F _(dn−1, 2n−3) ,

F_(dn, 2n−1) −F _(dn−1, 2n−2) ,

F_(dn, 2n])  (26)

(Step 5)

By representing the foregoing as follows, a Preisach distributionfunction table showing the hysteretic characteristic model is obtained:$\begin{matrix}{{\eta \left( {x_{i},x_{d}} \right)} = \begin{bmatrix}{\eta \left( {{- X_{n}},X_{n}} \right)} & {\eta \left( {{- X_{n - 1}},X_{n}} \right)} & \cdots & {\eta \left( {{- X_{1}},X_{n}} \right)} & {\eta \left( {X_{1},X_{n}} \right)} & {\eta \left( {X_{1},X_{n}} \right)} & \cdots & {\eta \left( {X_{n - 1},X_{n}} \right)} \\{\eta \left( {{- X_{n}},X_{n - 1}} \right)} & ⋰ & \vdots & \vdots & \vdots & \vdots & \ddots & 0 \\\vdots & \cdots & {\eta \left( {{- X_{2}},X_{2}} \right)} & {\eta \left( {{- X_{1}},X_{2}} \right)} & {\eta \left( {X_{0},X_{2}} \right)} & {\eta \left( {X_{1},X_{2}} \right)} & 0 & \quad \\{\eta \left( {{- X_{n}},X_{1}} \right)} & \cdots & {\eta \left( {{- X_{2}},X_{1}} \right)} & {\eta \left( {{- X_{1}},X_{1}} \right)} & {\eta \left( {X_{0},X_{1}} \right)} & 0 & \quad & \quad \\{\eta \left( {{- X_{n}},X_{0}} \right)} & \cdots & {\eta \left( {{- X_{2}},X_{0}} \right)} & {\eta \left( {{- X_{1}},X_{0}} \right)} & 0 & \quad & \quad & \vdots \\{\eta \left( {{- X_{n}},{- X_{1}}} \right)} & \cdots & {\eta \left( {{- X_{2}},{- X_{1}}} \right)} & 0 & \quad & ⋰ & \quad & \quad \\\vdots & \ddots & 0 & \quad & \quad & \quad & ⋰ & \quad \\{\eta \left( {{- X_{n}},{- X_{n - 1}}} \right)} & 0 & \quad & \quad & \cdots & \quad & \quad & 0\end{bmatrix}} & (27)\end{matrix}$

Each of FIGS. 12 and 13 shows a process when the amplitude ofdisplacement is increased from X₁ to X₂. This process will now bedescribed.

[1] If a load increase ΔF is set as F_(i11) when the displacement isincreased from −X₁ to 0, and loads are respectively F⁻¹² and F₁₁ at −X₁and 0, F_(i11) becomes a hatched part (1)+(3) of FIG. 12(a).

F _(i11) =F ₁₁ −F ⁻¹²  (28)

[2] If a load increase is F_(i12) when the displacement is increased toX₁, a load is F₁₂ at X₁, F_(i12) becomes (1)+(2)+(3), and the value of(2) is decided as shown in FIG. 12(b).

 F _(i12) =F ₁₂ −F ₁₁  (29)

η(X ₀ , X ₁)=F _(i12)  (30)

[3] Then, if a load decrease is F_(d11) when the displacement isdecreased from X₁ to 0, and a load is F⁻¹¹ at 0, F_(d11) becomes (1)+(2)of FIG. 12(c), and the value of η(−X₁, X₀) is decided.

F _(d11) =−(F ⁻¹¹ −F ₁₂)  (31)

η(−X ₁ , X ₀)=F _(d12)  (32)

Then, the value of η(−X₁, X₁) of (1) is obtained from the obtained η(X₀,X₁) of (2).

F _(d12) =−(F ⁻¹² −F ⁻¹¹)  (33)

η(−X ₁ , X ₁)=F _(d11) −F _(i12)  (34)

[4] Similarly, the following process is carried out for the amplitudeX₂(X₂>X₁) of a next one cycle.

As shown in FIG. 13(d), if loads are respectively F⁻²⁴ and F₂₁ at −X₂and −X₁, a load change F_(i21) when the displacement is increased from−X₂ to −X₁ becomes the added value of (4) to (7).

F _(i21) =F ₂₁ −F ⁻²⁴  (35)

[5] If a load is F₂₂ at the displacement of 0, a load change F_(i22)when the displacement is increased from −X₁ to 0 becomes the added valueof (1), (3) and (4) to (8). Since the added value of (1), (3), and (4)to (7) can be obtained, η(−X₂, −X₁) of (8) is obtained therefrom asshown in FIG. 13(e).

F _(i22) =F ₂₂ −F ₂₁  (36)

η(−X2, −X1)=F _(d22) −F _(i11)  (37)

[6] If a load is F₂₃ at X₁, a load change F_(i23) when the displacementis increased from 0 to X₁ becomes the added value of (1) to (9). Sincethe added value of (1) to (8) can be obtained, the value η(−X₂, X₀) of(9) is decided as shown in FIG. 13(f).

F _(i23) =F ₂₃ −F ₂₂  (38)

η(−X ₂ , X ₀)=F _(i23) −F _(i12)  (39)

[7] If a load is F₂₄ at X₂, a load change F_(i23) when the displacementis increased from X₁ to X₂ becomes the added value of (1) to (10). Sincethe added value of (1) to (9) can be obtained, the value η(−X₂, X₁) of(10) is obtained therefrom as shown in FIG. 13(g).

F _(i24) =F ₂₄ −F ₂₃  (40)

η(−X ₂ , X ₁)=F _(i24)  (41)

[8] If a load is F⁻²¹ at X₁, a load change F_(d21) when the displacementis decreased from X₂ to X₁ becomes the added value of (7) to (10). Sincethe added value of (8) to (10) can be obtained, the value η(−X₂, X₂) of(7) is obtained therefrom as shown in FIG. 13(h).

F _(d21)=−(F ⁻²¹ −F ₂₄)  (42)

η(−X ₂ , X ₂)=F _(d21)−(η₈+η₉+η₁₀)  (43)

[9] If a load is F⁻²² at the displacement of 0, a load change F_(d22)when the displacement is decreased from X1 to 0 becomes the added valueof (1), (2), (6) and (7) to (10). Since the added value of (1), (2) and(7) to (10) can be obtained, the value η(−X₁, X₂) of (6) is obtainedtherefrom as shown in FIG. 13(i).

F _(d22)=−(F ⁻²² −F ⁻²¹)  (44)

η(−X ₁ , X ₂)=F _(d22) −F _(d11)  (45)

[10] If a load is −F₂₃ at −X₁, a load change F_(d23) when thedisplacement is decreased from 0 to −X₁ becomes the added value of (1)to (3) and (5) to (10). Since the added value of (1) to (3) and (6) to(10) is obtained, the value η(X₀, X₂) of (5) can be obtained therefromas shown in FIG. 13(j).

F _(d23)=−(F ⁻²³ −F ⁻²²)  (46)

η(X ₀ , X ₂)=F _(d23) −F _(d12)  (47)

[11] If a load is −F₂₄ at −X₂, a load change F_(d24) when thedisplacement is decreased from −X₁ to −X₂ becomes the added value of (1)to (10). Since the added value of (1) to (3) and (5) to (10) can beobtained, the value η(X₁, X₂) of (4) is obtained therefrom as shown inFIG. 13(k).

F _(d24)=−(F ⁻²⁴ −F ⁻²³)  (48)

η(X ₁ , X ₂)=F _(d24)  (49)

In the above-described process, amplitude is increased as X₁, X₂, X₃, .. . , (X₁<X₂<X₃<. . . ) and, by repeating the processing until desireddisplacement amplitude is reached, a distribution function table can beformed.

When a hysteretic characteristic model is made based on the data of thehysteretic characteristic obtained from the vibration test by the methodof the invention, a more detailed hysteretic characteristic model can bemade by two-dimensionally interpolating the Preisach distributionfunction table using a spline function or the like.

By using the method of step 1, displacements x were provided from −1 to1 by 0.05 to the hysteretic characteristic of the Ramberg-Osgood typerepresented by (50), (51), and (52), α=4, β=1 and γ=2 were set, and adistribution function table was formed from loads F obtained therefrom.The result is shown in FIG. 14. In addition, the hystereticcharacteristic obtained from the distribution function is shown in FIG.15.

Skeleton curve: x=αF+βF ^(γ)  (50)

$\begin{matrix}{{\text{Force increase line:}\frac{x + x_{0}}{2}} = {{\alpha \frac{F + F_{0}}{2}} + {\beta \left( \frac{F + F_{0}}{2} \right)}^{\gamma}}} & (51)\end{matrix}$

$\begin{matrix}{{\text{Force reduction line:}\frac{x - x_{0}}{2}} = {{\alpha \frac{F - F_{0}}{2}} + {\beta \left( \frac{F - F_{0}}{2} \right)}^{\gamma}}} & (52)\end{matrix}$

Next, description will be made on a frequency response characteristicanalysis by an equivalent linearization method of a system having thehysteretic characteristic obtained from the above-described Preisachmodel. In this case, a solution is found to multi-degree-of-freedomnonlinear simultaneous equations by a procedure shown in FIG. 16.

(Step 6)

An equation of motion in the nonlinear system with hysteresis inrigidity is represented as follows.

m{umlaut over (x)}+f(x, {dot over (x)})=F cos ωt  (53)

Hereinafter, {umlaut over (x)} is denoted by Ax; and {dot over (x)} byBx.

Considering that the hysteretic characteristic is represented by f(x,Bx), a solution to the equation of motion is set approximately asfollows:

x=X cos φ, φ=ωt+θ  (54)

An output f(x, Bx) with respect to the hysteretic characteristic becomesa 2π-periodic function with respect to φ.

The output of the Preisach model with respect to an input represented bythe equation (54) is obtained by each of (13), (14), (15) and (16) (seeFIGS. 4(a) to 4(c)).

(Step 7)

The above output is subjected to Fourier development by the FFT, and aharmonic content left is represented as follows:

f(X cos φ, −Xω sin θ)=A ₁(X)cos φ+B ₁(X)sin φ  (55) $\begin{matrix}{{A_{1}(X)} = {\frac{1}{2\pi}{\int_{0}^{2\pi}{{f\left( {{X\quad \cos \quad \varphi},{{- X}\quad \omega \quad \sin \quad \varphi}} \right)}\cos \quad \varphi \quad {\varphi}}}}} & (56) \\{{B_{1}(X)} = {\frac{1}{2\pi}{\int_{0}^{2\pi}{{f\left( {{X\quad \cos \quad \varphi},{{- X}\quad \omega \quad \sin \quad \varphi}} \right)}\sin \quad \varphi \quad {\varphi}}}}} & (57)\end{matrix}$

If the equation of motion is replaced by an equivalent linear systemusing the above, it is represented as follows:

m{umlaut over (x)}+c _(e)(X){dot over (x)}+k _(e)(X)x=F cos ωt  (58)

Here, equivalent rigidity and an equivalent damping coefficient arerepresented as follows: $\begin{matrix}{{{k_{e}(X)} = \frac{A_{1}(X)}{X}},\quad {{c_{e}(X)} = {- \frac{B_{1}(X)}{X\quad \omega}}}} & (59)\end{matrix}$

By calculating the equation (59) and then solving the equation (58), anapproximate solution can be found.

Since the output (55) obtained from the hysteretic characteristicrepresented by the Preisach model with respect to the input of theequation (54) takes discrete values, Fourier coefficient (56), (57)thereof are numerically obtained. Thus, the equation (59) is obtained.

In the described embodiment, k_(e)(X) and c_(e)(X) in a range of Xavailable in the equation (59) are obtained beforehand. The range of X,and k_(e)(X) and c_(e)(X) are represented as follows:

X=[X ₁ , X ₂ , . . . , X _(n)]  (60)

k _(e)(X)=[k _(e)(X ₁), k _(e)(X ₂), . . . , k _(e)(X _(n))]  (61)

c _(e)(X)=[c _(e)(X ₁), c _(e)(X ₂), . . . , c _(e)(X _(n))]  (62)

Since it is not necessary to calculate the equation (55) for theequation (54) obtained at each convergent calculation, calculationefficiency is improved.

(Step 8)

From the equation (58), amplitude and phase of this harmonic vibrationare obtained from the following equation. $\begin{matrix}{{{X\sqrt{\left\{ {{k_{e}(X)} - {m\quad \omega^{2}}} \right\}^{2} + \left\{ {{c_{e}(X)}\omega} \right\}^{2}}} - F} = 0} & (63) \\{\varphi = {\tan^{- 1}\quad \frac{{c_{e}(X)}\omega}{{k_{e}(X)} - {m\quad \omega^{2}}}}} & (64)\end{matrix}$

This method is applied to the multi-degree-of-freedom system. If [M] isset as a mass matrix; [C_(e)(X)] as a matrix composed of equivalentdamping coefficients Ce(X); and [K_(e)(X)] as a matrix composed ofequivalent rigidity coefficients K_(e)(X), then the equation of motionis represented as follows:

└−ω² [M]+iω[C _(e)(X)]+[K _(e)(X)]┘X−F=[0]  (65)

X and F become vectors as follows:

X=(X ₁ , X ₂ , . . . , X _(n))  (66)

F=(F ₁ , F ₂ , . . . , F _(n))  (67)

Hereinafter, X is denoted by VX, and F by VF.

(Step 9)

When the Newton-Raphson method is applied to these simultaneousequations, C_(e)(X) and K_(e)(X) respectively constituting [C_(e)(X)]and [K_(e)(X)] become functions of displacement amplitude X.Accordingly, to obtain a Jacobian matrix thereof, it is necessary toobtain a derived function of the equivalent damping coefficient c_(e)(X)and the equivalent rigidity coefficient k_(e)(X).

Thus, according to the embodiment, interpolation is carried out based ona spline function of a piecewise polynomial for the obtained [C_(e)(X)]and [K_(e)(X)], and a derived function is obtained by differentiatingthis spline function.

If α_(i) is a constant coefficient, and p(x) is an mth order polynomial,an mth order spline function having N points of contact x=X₀, X₁, . . ., X_(N−1) is represented by the following equation (68). $\begin{matrix}{{s(x)} = {{p(x)} + {\sum\limits_{i = 0}^{N - 1}{\alpha_{i}\left( {x - X_{i}} \right)}_{+}^{m}}}} & (68)\end{matrix}$

Since the piecewise polynomial is differentiable, the derivative of themth order spline function is obtained by the following equation (SeeFIG. 17). $\begin{matrix}{{\frac{}{x}\left( {x - X_{i}} \right)_{+}^{m}} = {m\left( {x - X_{i}} \right)}_{+}^{m - 1}} & (69)\end{matrix}$

By using this method, the derivatives of the factors of [C_(e)(X)] and[K_(e)(X)] are obtained, a Jacobian matrix f′(VX) of the equation ofmotion f(VX) is formed. $\begin{matrix}{{f^{\prime}(X)}^{- 1} = \begin{bmatrix}\frac{\partial f_{1}}{\partial X_{1}} & \frac{\partial f_{1}}{\partial X_{2}} & \ldots & \frac{\partial f_{1}}{\partial X_{n}} \\\frac{\partial f_{2}}{\partial X_{1}} & \frac{\partial f_{2}}{\partial X_{2}} & \ldots & \frac{\partial f_{1}}{\partial X_{n}} \\\vdots & \vdots & ⋰ & \vdots \\\frac{\partial f_{n}}{\partial X_{1}} & \frac{\partial f_{n}}{\partial X_{2}} & \ldots & \frac{\partial f_{n}}{\partial X_{n}}\end{bmatrix}} & (70)\end{matrix}$

(Step 10)

If a Newton-Raphson method is applied to the equation (65), thefollowing is established:

└−ω² [M]+iω[C _(e)(X)]+[K _(e)(X)]┘X−F=f(X)=0  (71)

An approximate solution X^(k+1) of K+1 th time with respect to X^(k) ofKth time is represented as follows:

X ^(k+1) =X ^(k) −f′(X ^(k))⁻¹ /f(X ^(k)), k=0, 1, 2, . . .   (72)

A solution is found by repeating the equation (72). By carrying out thiscalculation for each ω of the equation (65), a frequency responsecharacteristic like that shown in FIG. 18 is obtained.

What is claimed is:
 1. A method for analyzing nonlinear restoring forcecharacteristic with hysteresis of a machine structure system,comprising: performing approximate mode analysis by representing saidnonlinear restoring force characteristic with hysteresis by anequivalent linear characteristic, when a nonlinear restoring forcecharacteristic of a machine structure system, having force-displacementhysteresis, is analyzed, wherein the representation of the nonlinearrestoring force characteristic by said equivalent linear characteristicincludes the steps of: first obtaining a basic nonlinear restoring forcecharacteristic experimentally or simulatively based on a relationbetween an input and an output of an actual system; then obtaining atransmission characteristic, equivalent rigidity, and equivalent viscousdamping by numerical analysis based on a relation between a sine waveinput with respect to said basic nonlinear restoring forcecharacteristic and a fundamental wave frequency component of an outputof the sine wave input distorted after passing through said basicnonlinear restoring force characteristic; and obtaining a frequencyresponse characteristic of the entire machine structure system byuniting said transmission characteristic, said equivalent rigidity andsaid equivalent viscous damping with a multi-degree-of-freedom equationof motion.
 2. A method for analyzing nonlinear restoring forcecharacteristic with hysteresis of a machine structure system,comprising: performing approximate mode analysis by representing saidnonlinear restoring force characteristic with hysteresis by anequivalent linear characteristic, when a nonlinear restoring forcecharacteristic of a machine structure system, having force-displacementhysteresis, is analyzed, wherein the representation of the nonlinearrestoring force characteristic by said equivalent linear characteristicincludes the steps of: first obtaining a basic nonlinear restoring forcecharacteristic experimentally or simulatively based on a relationbetween an input and an output of an actual system; then obtaining atransmission characteristic, equivalent rigidity, and equivalent viscousdamping by numerical analysis based on a relation between a sine waveinput with respect to said basic nonlinear restoring forcecharacteristic and a fundamental wave frequency component of an outputof the sine wave input distorted after passing through said basicnonlinear restoring force characteristic; obtaining a frequency responsecharacteristic of the entire machine structure system by uniting saidtransmission characteristic, said equivalent rigidity and saidequivalent viscous damping with a multi-degree-of-freedom equation ofmotion; and wherein when a solution is found to a nonlinear equation ofmotion of said machine structure system by Newton-Raphson method, datais obtained on equivalent rigidity of a nonlinear factor, and then thedata of the equivalent rigidity of the nonlinear factor is inerpolatedby a piecewise polynomial.